# -*- coding: utf-8 -*-

'''
Created on Jun 1, 2011

@author: Peter Harrington
'''
from numpy import *
import matplotlib
import matplotlib.pyplot as plt


# 读数据，默认tab分隔符
def loadDataSet(fileName, delim='\t'):
    # 打开文件
    fr = open(fileName)
    # 按行读取
    stringArr = [line.strip().split(delim) for line in fr.readlines()]
    # 转成float
    datArr = [map(float, line) for line in stringArr]
    return mat(datArr)


# *****PCA
# dataMat:样本数据
# topNfeat：最大特征值所取的个数
def pca(dataMat, topNfeat=9999999):
    # 求每个特征的均值（按列）
    # axis = 0：压缩行，对各列求均值，返回 1* n 矩阵
    # axis =1 ：压缩列，对各行求均值，返回 m *1 矩阵
    meanVals = mean(dataMat, axis=0)
    # 去均值，数据居中
    meanRemoved = dataMat - meanVals  # n*d - 1*d 结果是每一行都会减
    # 协方差矩阵
    covMat = cov(meanRemoved, rowvar=0)
    # 求特征值和特征向量
    eigVals, eigVects = linalg.eig(mat(covMat))
    # 特征值从小到大排序，返回位置索引index
    eigValInd = argsort(eigVals)  # 升序
    # 取最大的topNfeat个特征值的位置index
    eigValInd = eigValInd[:-(topNfeat + 1):-1]
    # 取特征值对应的特征向量
    redEigVects = eigVects[:, eigValInd]
    # 去均值数据 转化到新的低维空间
    lowDDataMat = meanRemoved * redEigVects
    # 为了使降维数据与原始数据在同一空间对比，需重构回原空间
    reconMat = (lowDDataMat * redEigVects.T) + meanVals
    return lowDDataMat, reconMat


def replaceNanWithMean():
    # 读数据
    datMat = loadDataSet('secom.data', ' ')
    # 特征数
    numFeat = shape(datMat)[1]
    # 遍历每个特征值
    for i in range(numFeat):
        # 每列求非NaN的均值
        meanVal = mean(datMat[nonzero(~isnan(datMat[:, i].A))[0], i])
        # 每列中值NaN替换为该列的均值
        datMat[nonzero(isnan(datMat[:, i].A))[0], i] = meanVal
    return datMat


if __name__ == '__main__':
    # #=================简单二维数据降维
    # topNfeat=1 #取特征值个数
    # #读取数据
    # dataMat=loadDataSet('testSet.txt')
    # #PCA降维
    # lowDMat,reconMat=pca(dataMat,topNfeat)
    #
    # #原始空间中绘图
    # fig = plt.figure()
    # ax = fig.add_subplot(111)
    # #显示原始数据
    # ax.scatter(dataMat[:, 0].tolist(), dataMat[:, 1].tolist(), marker='^', s=80)  # tolist()是后面加上去的，不加报错
    # #显示降维后数据
    # ax.scatter(reconMat[:,0].tolist(), reconMat[:,1].tolist(), marker='o', s=20, c='red')
    # plt.show()

    # ==================半导体数据的特征值与特征向量
    dataMat = replaceNanWithMean()
    meanVals = mean(dataMat, axis=0)
    meanRemoved = dataMat - meanVals  # 去均值
    covMat = cov(meanRemoved, rowvar=0)  # 协方差
    eigVals, eigVects = linalg.eig(mat(covMat))  # 特征值与特征向量
    # print eigVals
    #
    # ***选取的最大的几个特征值能量占总能量的比例
    eigValInd = argsort(eigVals)  # sort, sort goes smallest to largest
    # 特征值从大到小排序
    eigValInd = eigValInd[::-1]
    sortedEigVals = eigVals[eigValInd]
    total = sum(sortedEigVals)  # 所有特征值求和，即总方差
    varPercentage = sortedEigVals / total * 100  # 各个方差百分比
    print varPercentage[:20]

    fig = plt.figure()
    ax = fig.add_subplot(111)
    # 显示前20个特征值各自占的比重
    ax.plot(range(1, 21), varPercentage[:20], marker='^')
    plt.xlabel('Principal Component Number')
    plt.ylabel('Percentage of Variance')
    plt.show()